Many-body Anderson localization 
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We show, using quasi-exact numerical simulations, that Anderson localization of one-dimensional 
particles in a disordered potential survives in the presence of attractive interaction between par- 
ticles. The localization length of the composite particle can be computed analytically for weak 
disorder and is in good agreement with the quasi-exact numerical observations using Time Evolving 
Block Decimation. Our approach allows for simulation of the entire experiment including the final 
measurement of all atom positions. 
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Anderson localization (AL) has been widely investi- 
gated in the last 50 years 0, Q ■ The possibility of ob- 
serving directly localization of the wavefunction in cold 
atomic gases - in contrast with solid-state systems where 
transport properties such as the conductance are usu- 
ally measured - lead to a recent revival of the interest 
in localization properties in general, and in AL in par- 
ticular. AL is characterized by the inhibition of trans- 
port in a quantum system, whose classical counterpart 
behaves diffusively. It is accompanied by an exponen- 
tial localization of eigenstates in the configuration space, 
^(r)! 2 oc exp(— |r|/C), where ( is the localization length. 
AL is due to the interference between various multiple 
scattering paths which favors the return of the particle 
to its initial position and thus decreases its probability 
to travel a long distance. As the geometry of these paths 
depends on the system dimension, so AL features de- 
pend on the dimension too. For one-dimensional (ID) 
systems, AL is a generic single particle behavior even for 
very small disorder when the particle "flies" above the 
potential fluctuations. 

Being based on quantum destructive interference at 
long distance, AL is very sensitive to phenomena affecting 
the phase coherence of the wavefunction. For example, 
spontaneous emission from an excited state - or, in solid- 
state physics, phonon-electron scattering - will spoil the 
phase of the atomic (electronic in the solid-state context) 
wavefunction, leading to the destruction of AL and the 
restoration of a diffusive transport, albeit with a reduced 
diffusion constant at long times [!, Hj]. 

A fundamental question is to understand how inter- 
action between particles affects AL. If particles can be 
considered as uncorrelated, particle-particle interaction 
just acts like a decoherence process and will destroy AL. 
More interesting situations appear when the various par- 
ticles cannot be treated independently, for example be- 
cause they are prepared in a correlated state or because 
the inter-particle dynamics creates strong correlations. 



Presently, there is no consensus on the possible existence 
and properties of many-body localization. Some results 
suggest that AL survives for few-body systems, although 
with a modified localization length Q; studies of cold 
bosonic systems in the mean-field regime show that AL is 
destroyed and replaced by a sub-diffusive behavior 0, 0] ■ 
There are even predictions that AL survives at finite tem- 
perature in the thermodynamic limit, in the presence of 
interactions ||. In this Letter, we show, using a specific 
example, that ID AL survives in the presence of attrac- 
tive interactions and is even a rather robust phenomenon. 
The novelty of our approach is that it uses a quasi-exact 
numerical scheme to solve the full many-body problem in 
the presence of disorder. Here, quasi-exact means that 
all numerical errors can be controlled and reduced below 
an arbitrary value, just at the cost of increased compu- 
tational resources. The big advantage of this approach is 
not to rely on neglecting a priori some physical processes, 
which is very dangerous in this many-body quantum con- 
text where any small uncontrolled perturbation may lead 
to decoherence and destruction of AL. 

Atomic matter waves have several advantages that 
made possible an unambiguous demonstration of single 
particle AL in ID 0, |T3| : atom-atom interaction can be 
reduced either by diluting the atomic gas or using Fes- 
hbach resonances, ensuring a very long coherence time 
of the atomic matter wave; the spatial and temporal or- 
ders of magnitudes are very convenient allowing a direct 
spatio-temporal visualization of the dynamics; all micro- 
scopic ingredients are well controlled; a disordered po- 
tential can be created by using the effective potential 
induced by a far detuned optical speckle. Here we shall 
not avoid interactions, on the contrary, interactions form 
a central issue of the study. 

We consider a system of N identical bosonic atoms in a 
ID system, in the regime of attractive two-body interac- 
tions. We assume the dilute regime where the atom-atom 
average distance is larger than the scattering length and 
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take the low-energy limit where the interaction can be 
modeled by a (negative) (5-potential. 

The many-body Hamiltonian can be written, using the 
standard second quantization formalism (assuming unit 
mass for the particles and taking H = 1): 



H = J dz ft(z) 
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available computer resources, a maximum of N =25 par- 
ticles could be included in the calculations; similar results 
are obtained for 10 particles. We use the soliton size £ as 
the unit of length; consequently the time unit is £ 2 . We 
choose the harmonic trap shallow enough not to distort 
the soliton shape, but still strong enough to confine its 
center of mass over a distance only slightly larger than its 

size, that is Aq = y§£- This relates to choosing the fre- 

(1) quency ui of the trapping harmonic potential ui 2 z 2 /2 such 



where g < is the strength of the atom-atom interaction 
and V(z) the external potential. 

In the absence of disordered potential, the ground state 
of this system is known exactly, thanks to the Bethe 
ansatz [11] . For large N, it is conveniently described as a 
bright soliton, a composite particle whose degrees can be 
separated between internal degrees of freedom and two 
external degrees of freedom: a global phase which is 
irrelevant as long as a single soliton is considered, and 
the position q of the center of mass of the soliton. The 
particle density of the soliton - normalized to the total 
number N of particles - is given by \4>o( z ~ q)\ 2 where <fio 
is the ground state of the mean-field Hamiltonian: 



<Po{z - q) = 



2£ cosh [(z -</)/£]" 



(2) 



£ = —-fig is the characteristic size of the soliton. The 
associated chemical potential is [i = —N 2 g 2 /8. In this 
mean-field approach, the center of mass q is a simple 
parameter. For the full many-body system, q must be 
thought of as the quantum position operator of the soli- 
ton, the composite particle formed by the N particles. 

With an additional disordered potential V(z) the 
many-body problem cannot be solved exactly. One must 
rely on quasi-exact numerical approaches. A convenient 
way is to discretize the continuous Hamiltonian, eq. (fl). 
over a discrete lattice ■ The resulting Hamiltonian has 
a standard Bose-Hubbard form [l3[ enabling us to use 
well developed tools for ID lattice systems. The ground 
state as well as the dynamics may be quasi-exactly stud- 
ied using the Time Evolving Block Decimation (TEBD) 
algorithm [Til . HEl ]. essentially equivalent to the time- 
dependent Density Matrix Renormalization Group ap- 
proach (lil . 03] ■ The key point is to describe at any time 
the state of the system as a so-called Matrix Product 
State (MPS) (see |llE3) ■ 

The numerical simulation proceeds as follows: in a first 
step, we compute the ground state of N interacting par- 
ticles in the presence of a harmonic trap (to have the ini- 
tial soliton localized at the center) using imaginary time 
propagation with TEBD. In a second step, we use this 
many-body state as the initial state for real time propa- 
gation with TEBD, in the absence of the harmonic trap, 
but in the presence of a weak speckle potential. With our 
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0.025/£ . When the harmonic 



trap is switched off and the laser speckle switched on, we 
expect the many-body composite particle (the soliton) 
to diffuse and eventually localize thanks to many-body 
AL. In order for the localization length to be reasonably 
short, we chose the strength of the external speckle po- 
tential comparable to the initial average energy of the 
soliton (w/4), that is Vq — 2.5 x 10~ 4 . We also chose the 
correlation length of the speckle potential <tq = 0.4£ to 
be significantly shorter than the soliton size, so that the 
effective one-body potential (see below) is generic and 
free of the peculiarities of the speckle potential (ljj ). 

Several sources of errors exist in the TEBD algorithm 
and must be controlled. The first one is due to the spa- 
tial discretization. Getting accurate results requires the 
discretization unit, 8, to be much smaller than both the 
soliton size £ and the de Broglie wavelength of the soli- 
ton 2n/k, with k of the order of the typical wavevector 
contained in the initial wavepacket. We use 8 = £/5; a 
twice smaller step produces slightly different quantitative 
results, but the difference is practically invisible on the 
scale of the figures shown below. The second source of 
error is temporal discretization of the evolution operator. 
We use the standard Trotter expansion [l3|, whose error 
can be controlled by varying the time step (St = 0.008£ 2 
is used). The third source of error is the truncation of 
the MPS at each step. This error is monitored through 
the so-called "discarded weight", that is the weight of the 
components which have to be discarded from the time 
evolved many-body state to keep it in MPS form with a 
fixed parameter \ ~ f ne number of bonds between sites 
fl3l | . This can be a serious problem when the entropy 
of entanglement grows as a function of time. It may 
even prevent calculation to be ran beyond some rather 
short time, especially when th e sy stem is significantly 
excited above the ground state [l9|. In our case, the en- 
ergy of the many-body wavepacket is small and we have 
not observed any significant growth of the entropy of en- 
tanglement when AL sets in. This made it possible to 
run the full calculation up to time t = 4000. That has 
been tested using a larger variational set representing the 



wavefunction [13 1 



Figure [T] shows the particle density in configuration 
space obtained at three different times t = 1000, t = 2000 
and the final time t = 4000, averaged over 96 realizations 
of the disorder. In the absence of disordered potential, 



FIG. 1: (Color on line) Spatial atomic density of an initially 
localized bright soliton (dotted line), under the influence of 
a disordered speckle potential, using the quasi-exact many- 
body TEBD algorithm: t=1000 (dashed line), 2000 (dash- 
dotted line), 4000 (solid line). After an initial expansion, the 
atomic density freezes at long times, a signature of many-body 
Anderson localization. Parameters are given in the text. 



the initial wavepacket is expected to spread. Single par- 
ticle physics gives the characteristic time for this spread- 
ing, t = 1/lo = 40, much shorter than the time scale in 
the figure. At £=1000, one clearly sees that the central 
part of the wavepacket is already more or less localized 
while the ballistic front for \z\ > 40 has not yet been 
scattered and keeps a Gaussian shape. This corresponds 
to the wavepacket components with the highest energy 
and consequently the longest localization length. AL has 
already setup at £=2000 and does no longer evolve fur- 
ther, compare with i=4000. Therefore, Fig. Q] provides 
an evidence for many-body AL taking place in a quasi- 
exact full many-body numerical simulation. At the final 
time (100 times the characteristic spreading time), we do 
not observe any indication that AL could be destroyed. 

The description of the final state as a MPS makes it 
possible to compute easily more complicated quantities, 
such as correlation functions. The simplest one is the 
one-body density matrix (tfj(z)ip^ (z')} , shown in Fig. [2] 
It clearly displays extremely strong correlations between 
positions z and z', reinforcing the observation that AL 
probably survives far beyond t =4000. The interpreta- 
tion is simple in terms of bright soliton: all atoms are 
grouped in a soliton of size £ = 1 , but the center of mass 
of the soliton itself is widely spread. This has an impor- 
tant consequence: the largest eigenvalue of the one-body 
density matrix - a value often used as a quantitative cri- 
terion for Bose- Einstein condensation [20| - is here 0.14, 
much smaller than unity; in contrast, the value at t = 
is 0.84. Thus, while the initial state can be considered 
as a true condensate, the temporal dynamics destroys 
condensation; any description of our many-body system 
using the mean field theory via the Gross-Pitaevskii equa- 
tion (which by construction describes a 100% condensate) 
must fail. In other words, our many-body AL is neces- 
sarily beyond the mean field description. 

The structure of the MPS is simple enough to make 
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FIG. 2: (Color on line) One-body density matrix of the many- 
body state at t = 4000 (see fig. [1]) in configuration space. It 
is strongly concentrated along the diagonal z — z' with a 
transverse width of the order of the soliton size £. This is a 
direct proof that the many-body system can be described by 
a compact composite particle, a bright soliton, whose center 
of mass is widely spread. 



more complicated properties easily measurable. For ex- 
ample, it is possible to simulate the measurement of the 
particle positions. Indeed, the position of the first "mea- 
sured'' particle can be chosen randomly following the 
atomic density. Once the position z\ is chosen, it is easy 
to project the MPS state using the annihilation 
operator on another MPS state containing only N — 1 
particles. Iteration of the process allows us to measure 
all particle positions, see inset of Fig. [3] and [T^j . In this 
way, we extract both the position of the center of mass 
of the soliton and the atomic density relative to the cen- 
ter of mass. The later quantity is shown in figure [3] for 
time t = 4000 in comparison with the analytic predic- 
tion, eq. @. The agreement is excellent, showing that 
the internal structure of the soliton is fully preserved for 
a long time, even after AL has completely set in. The 
small difference is a 1/N finite size effect. 

It is possible to write an approximate effective Hamil- 
tonian describing the dynamics of the center of mass of 
the soliton |2ll423| . The internal degrees of freedom of 
the bright soliton are gapped, with an energy gap equal 
to —(j, = N 2 g 2 /8, so that a weak external perturbation 
cannot populate internal excited states of the bright soli- 
ton; this strongly contrasts the situation of repulsive in- 
teraction where one can form dark solitons which are 
not gapped and consequently more sensitive to external 
perturbations [24[. Assuming a fixed soliton shape, the 
center of mass q of the N particles can be described as a 
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FIG. 3: (Color on line) Atomic density measured with respect 
to the center of mass of the N particles, at time t = 4000, after 
Anderson localization has set in. A large number of individ- 
ual "measurements" of the particle positions - two examples 
are shown in the inset - are performed for various disorder re- 
alizations, and the resulting average is displayed in the main 
panel. The agreement between the full many-body calcula- 
tion (solid line) and the prediction of eq. ((2| for the soliton 
shape (dashed line) is excellent. 



collective variable, leading to the effective Hamiltonian: 



q ~ 2N 



+ / dz \Mz-q)\ 2 v(z) 



(3) 



where p q is the momentum conjugate to q. This Hamilto- 
nian describes a composite particle with mass N evolving 
in a potential V(q) that is the convolution of the bare po- 
tential with the soliton envelope. 

This effective one-body theory is able to quantitatively 
predict the long time limit for the spatial density prob- 
ability of the soliton center of mass, see the detailed 
derivation and calculations in [2l|. Initially, the soliton 
is in the ground state of the harmonic trap (a Gaus- 
sian wavepacket) that, after the trap is switched off, ex- 
pands over a range of energies, each energy component 
being characterized by its own localization length. As a 
conse quen ce, although each component displays approx- 
imate [25l | exponential localization in the long time limit, 
the superposition displays approximate algebraic local- 
ization, as discussed in [2l| . 

In figure [4l we show comparison between the full many- 
body calculation and the effective one-body theory. For 
the latter theory, we numerically propagated the initial 
Gaussian wavepacket in presence of the effective potential 
V(q). Note that the many-body result is here plotted for 
the center of mass position, which can slightly differ from 
the atomic density, the latter being the convolution of the 
former by the soliton shape. At the scale of the figure, 
the soliton is extremely narrow so that both quantities 
are almost equal, compare with fig. [1] The agreement 
between the many-body and the one-body calculations is 
clearly excellent, validating the one-body approach. The 
many-body calculation produces a slightly larger density 
in the tail. This is due to the slightly too small size, 
[— 120£, 120£] used in the many-body calculation and re- 
flection of the head of the wavepacket on the boundaries. 




0.001. 



60 -40 -20 20 40 60 



Wq)\ - 
0.01- 



0.001 



-r. — i — i — i — i i r 



l 1 1 — r 




_i i i i i i i 



FIG. 4: (Color on line) Probability density for the center of 
mass of the soliton at time t — 4000, numerically computed 
using the full many-body quasi-exact dynamics (solid line, 
averaged over 96 disorder realizations) and the effective one- 
body theory (dashed line, averaged over 10000 realizations). 
The good agreement shows that the many-body problem ac- 
tually displays Anderson localization, which is correctly pre- 
dicted by a simple effective one-body theory for a composite 
particle. The dotted line in the lower panel (double logarith- 
mic plot) is the l/\q\ approximate prediction of the one-body 
effective theory. 



In fig. |4j we also show the l/\q\ leading behaviour pre 



dieted by the one-body theory, eq. (13) of Ref. [21|. It 
correctly predicts the observed behaviour. It does not 
aim at being quantitative for several reasons: the main 
ones are the existence of a subleading logarithmic term 
and the fact that this formula assumes a weak disorder, 
an assumption not fully satisfied here. 

To summarize, we have shown the existence of many- 
body Anderson localization for attractive bosons in the 
presence of a disordered potential. The claim is based 
on quasi-exact many-body numerical simulations using 
the TEBD algorithm, which incorporate all complicated 
phenomena that could spoil the internal phase coher- 
ence of the many-body composite object, a bright soli- 
ton, displaying Anderson localization. Moreover, we ob- 
tain excellent agreement between the many-body calcula- 
tion and a one-body effective theory, which goes beyond 
standard mean field theories such as the Gross-Pitaevskii 
equation. Our quasi-exact many-body approach allows 
for simulation of the entire experiment starting from the 
initial state in a harmonic strap till the destructive mea- 
surement of all atom positions. 
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The discretization of space to a chain of sites located at 
equally spaced positions q\ = 16 together with a 3-point 
discretization of the Laplace operator allows us to write 
the Hamiltonian in a tight-binding Bose-Hubbard form 

H = -J{a\ai + i + h.c.) + ^a\a\aiai + Vi ajaj 

(4) 

with J = 2^7, U = j and Vi = V(qi) (an additional 
trivial constant term 1/5 2 has been dropped). 

For numerical purposes the infinite space is restricted 
to the [-K5, K5] interval leading to a ID chain of M = 
2K + 1 discrete sites. N identical bosons are distributed 
on these M sites. A basis of the total Hilbert space can 
be built using the direct product of Fock states on each 
site \i\,i2 ■ ■ ■ ijvf ) with the constraint that the sum of oc- 
cupation numbers 53;=i m *J ^ s equal to N. The size of 
this basis increases exponentially with the system size, 
making its use unpractical for large many-body problems. 
Instead, we use a variational set of Matrix Product States 
(MPSs). An MPS is a state which can be written as: 



|V) = 



■nWM \ [1] r [2],i 2 r [M],i M | • 
"ai 1 «ioj ' • -1 a B _il I'll 



,lM 



) (5) 



where r^' 1 ' are matrices and A^ vectors. To describe 
exactly a generic state in terms of MPS, a large number 
(exponentially increasing with M) of ai values is needed. 
However, typical low energy states are only slightly en- 
tangled so that A^' =1 2 are rapidly decaying numbers, 
which allows for introduction of a cutoff \ m the sum 
over Greek indices above, resulting in tractable numeri- 
cal computations [ljj ■ For a ground state protected by a 
gap, the area theorem (2(| ensures that an efficient MPS 
representation exists. In the physical situation discussed 
by us, the area law has no direct applicability. 

The i\,..iM indices are in principle restricted to the 
[0, N] interval. In practice, since it is highly unlikely 
that all the bosons occupy a single site of the system, we 
lower the cutoff in the sums assuming some iV max < N. 
While the maximum average occupation number (at sites 
near the center of mass of the soliton) is N5/2t;=2.5, we 
found surprisingly that a high -/V max is needed for conver- 
gent results, while the x value may be kept relatively low. 
The converged calculations reported use 7V max =14, x=30 
yielding the internal Hilbert space dimension per site be- 
ing 450. The results has been compared with those with 
lower iV m ax as well as those for x=40 (for shorter times) 
to check that the results presented are fully converged. 

Referring the reader to original papers and reviews 

0, 

15 . ItJ for the details, let us mention only that the TEBD 
algorithm describes how the and A^ evolve in time 
under the influence of an Hamiltonian containing simple 



6 



terms local on each site as well as hoping terms of the 
type aia\ +l which transfer one particle from site I to site 
1 + 1. 

Being able to write the many-body state as a MPS has 
considerable advantages, especially if it is - like in our 
calculations - in the so-called canonical form [17j , For 
example, expectation values of local operators such as 
ajai or a\ai+\ involves only simple contractions on the 
local r"l tensors and A^ vectors. It makes it also possible 
to mimic the measurement process of particle positions 
as follows. The reduced density matrix on site I is 
easily constructed by contracting the tensor with the 
neighboring A'' -1 ' and A^ vectors: 

p [1] -= V W [l - 1] } 2 r \\W} 2 (6) 

ai-i,ai 

We then chose randomly the number of particles "mea- 
sured" on site I following the statistical populations, di- 
agonal elements of the on-site reduced density matrix. 
Once a given occupation number i is chosen, we project 



the MPS state onto the subspace with ii = i and nor- 
malize it. This involves only simple contractions on the 
local T tensors and A vectors, producing another MPS. 
The process can be iterated on all sites, and is partic- 
ularly simple if sites are scanned consecutively starting 
from one edge and propagating toward the other edge. 
It is simple to prove that the probability distribution of 
the measurements is independent of the order used for 
scanning the various sites. 

An individual "measurement" produces a single set of 
occupation numbers (ii, 12---im) (whose sum is of course 
TV) whose probability is exactly i2---iM\4>)\ 2 ■ By per- 
forming a series of "measurements", we can sample in- 
teresting physical quantities, such as the position of the 
soliton center of mass, (q) = J2i ISii/N, and the particle 
density with respect to this center of mass. Note that 
these quantities are hard to measure by other means as 
they involve correlation functions of high order (typically 
up to N) |22ll23. 



